#ifndef __lagrange_1d_hh__
#define __lagrange_1d_hh__

#include "line.hh"
#include "utils.hh"
#include "femlib_exceptions.hh"
// #include "unroll.hh"

#include "polylib_templated.hh"

#if defined (FADBADPP_GRADIENTS) || defined (FADBADPP_HESSIANS)
#include "fadiff.hh"
#include "badiff.hh"
#include "tadiff.hh"
#endif


// ALL Lagrange type 1-D basis functions: nth Order elements
template <unsigned int ORDER, typename T=double>
class LagrangeSegmentFiniteElement : public SegmentFiniteElement<ORDER, fastfem::Lagrange, T> 
{
public:

  typedef FEMTraits<LagrangeSegmentFiniteElement<ORDER, T> > Traits;

  LagrangeSegmentFiniteElement() : SegmentFiniteElement<ORDER, fastfem::Lagrange, T> ()
  {
    compute_vertices(vertices) ;
  }

  virtual ~LagrangeSegmentFiniteElement()
  {
  }

  /** virtual function which evaluates the values of all
       shape functions at a given point ip */
  virtual T evaluate_basis ( unsigned int , const Point<SPATIAL_DIMENSION, T>& ) ;

  /** virtual function which evaluates the values of all
       shape functions at all point ips */
  virtual void evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<T>& basis_vals ) ;

  /** virtual function which evaluates the values of all
       shape functions at a given point ip */
  template <typename U>
  static U evaluate_basis ( unsigned int , const Point<SPATIAL_DIMENSION, U>& ) ;

  /** virtual function which evaluates the values of all
       shape functions at all point ips */
  template <typename U>
  static void evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, U>& ips, std::vector<U>& basis_vals ) ;

  // void evaluate_basis_all_unrolled ( const PointList<SPATIAL_DIMENSION>& ips, std::vector<T>& basis_vals ) ;
  
  virtual void evaluate_gradient ( unsigned int , const Point<SPATIAL_DIMENSION, T>&, Point<SPATIAL_DIMENSION, T>& pt) ;

  virtual void evaluate_gradient_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<Point<SPATIAL_DIMENSION, T> >& grad_vals ) ;

  virtual void evaluate_hessian ( unsigned int , const Point<SPATIAL_DIMENSION, T>&, PointList<SPATIAL_DIMENSION, T>& pt) ;

  virtual void evaluate_hessian_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<PointList<SPATIAL_DIMENSION, T> >& grad_vals ) ;

  void reference_element_coords(std::vector<Point<SPATIAL_DIMENSION, T> >& nodes) ;  

  virtual int test() ;

#ifdef FADBADPP_GRADIENTS
  typedef fadbad::F<T> FADGradient;
  enum { FADGradient_Type = fastfem::FORWARD } ;
#endif

#ifdef FADBADPP_HESSIANS
  typedef fadbad::B<fadbad::F<T> > FADHessian;
  enum { FADHessian_Type = fastfem::BACKWARD_FORWARD } ;
#endif

  static const Point<SPATIAL_DIMENSION, T>* get_vertices() ;
  
protected:
  static Point<SPATIAL_DIMENSION,T> vertices[Traits::N_Nodes];
  static T weights[Traits::N_Nodes];
  static bool compute_vertices_called;
  
  static void compute_vertices(Point<SPATIAL_DIMENSION, T>[Traits::N_Nodes]) ;
};


template <unsigned int ORDER, typename T>
struct FEMTraits<LagrangeSegmentFiniteElement<ORDER, T> >
{
  
  enum { FE_Dimension = 1 } ;
  enum { FE_Family = fastfem::Lagrange } ;
  enum { FE_Order = ORDER } ;
  enum { N_Nodes = ORDER+1 } ;

  typedef T BaseType;
  typedef LagrangeSegmentFiniteElement<ORDER, T> FEType;
  typedef NodalFiniteElement<1, fastfem::Segment, ORDER, fastfem::Lagrange> ParentFEType;

};


// All inline functions

/** virtual function which evaluates the values of all
      shape functions at a given point ip */
template <unsigned int ORDER, typename T>
inline
T LagrangeSegmentFiniteElement<ORDER, T>::evaluate_basis (  unsigned int ib, const Point<SPATIAL_DIMENSION, T>& px)
{
  femlib_assert(ib < ORDER+1, "Invalid basis index request");

  const Point<SPATIAL_DIMENSION, T>* vertices = LagrangeSegmentFiniteElement<ORDER, T>::get_vertices() ;

  T bvalue = 1.0 ;
  for(unsigned int j=0; j < Traits::N_Nodes; ++j)
  {
    if (j != ib)
      bvalue = bvalue * (px.x() - vertices[j].x()) / (vertices[ib].x() - vertices[j].x());
  }

  return bvalue ;
}

/** virtual function which evaluates the values of all
      shape functions at all point ips */
template<unsigned int ORDER, typename T>
inline
void LagrangeSegmentFiniteElement<ORDER, T>::evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<T>& basis_vals )
{
  basis_vals.resize((ORDER+1)*ips.n_points, 1.0);

  const Point<SPATIAL_DIMENSION, T>* vertices = LagrangeSegmentFiniteElement<ORDER, T>::get_vertices() ;

  for (unsigned int k = 0, offset=0; k < ips.n_points; k++, offset+=(ORDER+1) )
  {
    const double& x = ips.x(k) ;
    for (unsigned int i = 0; i < Traits::N_Nodes; i++ )
      for (unsigned int j = 0; j < Traits::N_Nodes; j++ )
      {
        if ( j != i )
          basis_vals[i+offset] = basis_vals[i+offset] * ( x - vertices[j].x() ) / ( vertices[i].x() - vertices[j].x() );
      }
  }
  
  return ;
}


/*
template<unsigned int ORDER>
inline
void LagrangeSegmentFiniteElement<ORDER>::evaluate_basis_all_unrolled ( const PointList<SPATIAL_DIMENSION>& ips, std::vector<double>& basis_vals )
{
  int sum = 0;
  int tmp;
  //double *vert = this->vertices;
  
  fastfem::UnrollerL<0, Traits::N_Nodes>::step( [&] (int i) {
      fastfem::UnrollerL<0, Traits::N_Nodes>::step( [&] (int j) {
        if ( j != i )
          for (unsigned int k = 0; k < ips.n_points; k++ )
            basis_vals[i+k*(Traits::N_Nodes)] = basis_vals[i+k*(Traits::N_Nodes)] * ( ips.x(k) - this->vertices[j].x() ) / ( this->vertices[i].x() - this->vertices[j].x() );
        }
      )
    }
  );
}
*/

template <unsigned int ORDER, typename T>
inline
void LagrangeSegmentFiniteElement<ORDER, T>::evaluate_gradient (  unsigned int ib, 
                                                                const Point<SPATIAL_DIMENSION, T>& px, 
                                                                Point<SPATIAL_DIMENSION, T>& grad)
{
  femlib_assert(ib < Traits::N_Nodes, "Invalid basis index request");
  
#ifdef FADBADPP_GRADIENTS
  
  Point<SPATIAL_DIMENSION, FADGradient> pp(px);

  // Initialize variables x, y and z -> primary variables as per dimension
	FADGradient& x = pp[0];                 
#if (SPATIAL_DIMENSION > 1)
	FADGradient& y = pp[1];                 // Initialize variable x
#if (SPATIAL_DIMENSION > 2)
  FADGradient& z = pp[2];                 // Initialize variable x
#endif
#endif

  // Differentiate with respect to each primary variable
	x.diff(0,3);         // (index 0 of 3)

#if (SPATIAL_DIMENSION > 1)
	y.diff(1,3);         // (index 1 of 3)
#if (SPATIAL_DIMENSION > 2)
  z.diff(2,3);         // (index 2 of 3)
#endif
#endif
	
  // Evaluate function and derivatives
	FADGradient f = LagrangeSegmentFiniteElement<ORDER, T>::evaluate_basis<FADGradient>(ib, pp);
  // f.diff(0, 1) ;
  
  // Value of basis function at the requested point
	//double fval=f.x();
  
  // store the value of derivatives: df/dx into input variable 'grad'
	grad[0] = f.d(0);
#if (SPATIAL_DIMENSION > 1)
	grad[1] = f.d(1);
#if (SPATIAL_DIMENSION > 2)
  grad[2] = f.d(2);
#endif
#endif

#else
  
  const Point<SPATIAL_DIMENSION, T>* vertices = LagrangeSegmentFiniteElement<ORDER, T>::get_vertices() ;

  grad[0] = 0.0 ;
  for(unsigned int j=0; j < Traits::N_Nodes; ++j)
  {
    if (j != ib)
    {
      double prod = 1.0 ;
      for(unsigned int k=0; k < Traits::N_Nodes; ++k)
      {
        if (k != ib && k != j)
          prod = prod * (px.x() - vertices[k].x()) / (vertices[ib].x() - vertices[k].x());
      }
      grad[0] = grad[0] + prod / (vertices[ib].x() - vertices[j].x());
    }
  }
  
#endif

  return ;
}

template <unsigned int ORDER, typename T>
inline
void LagrangeSegmentFiniteElement<ORDER, T>::evaluate_gradient_all ( const PointList<SPATIAL_DIMENSION, T>& ips, 
                                                                      std::vector<Point<SPATIAL_DIMENSION, T> >& grad_vals )
{
  grad_vals.resize((Traits::N_Nodes)*ips.n_points);

  const Point<SPATIAL_DIMENSION, T>* vertices = LagrangeSegmentFiniteElement<ORDER, T>::get_vertices() ;

  for (unsigned int l = 0; l < ips.n_points; l++ )
  {
    for (unsigned int ib = 0; ib < Traits::N_Nodes; ib++ )
    {
      grad_vals[l*(Traits::N_Nodes)+ib][0] = 0.0 ;
      for(unsigned int j=0; j < Traits::N_Nodes; ++j)
      {
        if (j != ib)
        {
          T prod = 1.0 ;
          for(unsigned int k=0; k < Traits::N_Nodes; ++k)
          {
            if (k != ib && k != j)
              prod = prod * (ips.x(l) - vertices[k].x()) / (vertices[ib].x() - vertices[k].x());
          }
          grad_vals[l*(Traits::N_Nodes)+ib][0]= grad_vals[l*(Traits::N_Nodes)+ib][0] + prod / (vertices[ib].x() - vertices[j].x());
        }
      }
    }
  }
  
  return ;
}


template <unsigned int ORDER, typename T>
inline
void LagrangeSegmentFiniteElement<ORDER, T>::evaluate_hessian ( unsigned int ib, 
                                                              const Point<SPATIAL_DIMENSION, T>& px, 
                                                              PointList<SPATIAL_DIMENSION, T>& hess)
{
  femlib_assert(ib < Traits::N_Nodes, "Invalid basis index request");
  
#ifdef FADBADPP_HESSIANS
  
  Point<SPATIAL_DIMENSION, FADHessian> pp(px);

  // Initialize variables x, y and z -> primary variables as per dimension
	FADHessian& x = pp[0];
#if (SPATIAL_DIMENSION > 1)
	FADHessian& y = pp[1];                 // Initialize variable x
#if (SPATIAL_DIMENSION > 2)
  FADHessian& z = pp[2];                 // Initialize variable x
#endif
#endif

  // Differentiate with respect to each primary variable
	x.x().diff(0,3);         // (index 0 of 3)
#if (SPATIAL_DIMENSION > 1)
	y.x().diff(1,3);         // (index 1 of 3)
#if (SPATIAL_DIMENSION > 2)
  z.x().diff(2,3);         // (index 2 of 3)
#endif
#endif
	
  // Evaluate function and derivatives
	FADHessian f = LagrangeSegmentFiniteElement<ORDER, T>::evaluate_basis<FADHessian>(ib, pp);
  f.diff(0, 1) ;
  
  // Value of basis function at the requested point
	// double fval=f.x().x();
  // double dfdx=x.d(0).x();    // Value of df/dx
	// double dfdy=y.d(0).x();    // Value of df/dy
	// double dfdxdx=x.d(0).d(0); // Value of df/dx
	// double dfdxdy=x.d(0).d(1); // Value of df/dxdy
	// double dfdydx=y.d(0).d(0); // Value of df/dydx
	// double dfdydy=y.d(0).d(1); // Value of df/dy
  
  //hess.
  
  // store the value of derivatives: df/dx into input variable 'grad'
  for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
  {
    T* pval = hess.get(id);
    
    for (unsigned int jd=0; jd < SPATIAL_DIMENSION; ++jd)
      pval[jd] = pp[id].d(id).d(jd);
  }

#else

  // Reference: http://www.phys.ufl.edu/~coldwell/wsteve/FDERS/The%20Lagrange%20Polynomial.htm
  // m -> ib
  // i -> m
  // n -> j
  // j -> k

  const Point<SPATIAL_DIMENSION, T>* vertices = LagrangeSegmentFiniteElement<ORDER, T>::get_vertices() ;

  hess[0] = 0.0 ;
  for(unsigned int m=0; m < Traits::N_Nodes; ++m)
  {
    if (m != ib)
    {
      T dprod = 0.0 ;
      for(unsigned int j=0; j < Traits::N_Nodes; ++j)
      {
        if (j != m && j != ib)
        {
          T prod = 1.0 ;
          for(unsigned int k=0; k < Traits::N_Nodes; ++k)
          {
            if (k != m && k != j && k != ib)
              prod = prod * (px.x() - vertices[k].x()) / (vertices[ib].x() - vertices[k].x());
          }
          dprod = dprod + prod / (vertices[ib].x() - vertices[j].x());
        }
      }
      hess[0] = hess[0] + dprod / (vertices[ib].x() - vertices[m].x());
    }
  }

#endif

  return ;
}


template <unsigned int ORDER, typename T>
inline
void LagrangeSegmentFiniteElement<ORDER, T>::evaluate_hessian_all ( const PointList<SPATIAL_DIMENSION, T>& ips, 
                                                                  std::vector<PointList<SPATIAL_DIMENSION, T> >& hess_vals )
{
  hess_vals.resize((Traits::N_Nodes)*ips.n_points);

  const Point<SPATIAL_DIMENSION, T>* vertices = LagrangeSegmentFiniteElement<ORDER, T>::get_vertices() ;

  // Reference: http://www.phys.ufl.edu/~coldwell/wsteve/FDERS/The%20Lagrange%20Polynomial.htm
  // m -> ib
  // i -> m
  // n -> j
  // j -> k
  for (unsigned int l = 0; l < ips.n_points; l++ )
  {
    for (unsigned int ib = 0; ib < Traits::N_Nodes; ib++ )
    {
      T* pval = hess_vals[l*(Traits::N_Nodes)+ib].get(0);
      pval[0] = 0.0 ;
      for(unsigned int m=0; m < Traits::N_Nodes; ++m)
      {
        if (m != ib)
        {
          T dprod = 0.0 ;
          for(unsigned int j=0; j < Traits::N_Nodes; ++j)
          {
            if (j != m && j != ib)
            {
              T prod = 1.0 ;
              for(unsigned int k=0; k < Traits::N_Nodes; ++k)
              {
                if (k != m && k != j && k != ib)
                  prod = prod * (ips.x(l) - vertices[k].x()) / (vertices[ib].x() - vertices[k].x());
              }
              dprod = dprod + prod / (vertices[ib].x() - vertices[j].x());
            }
          }
          pval[0] = pval[0] + dprod / (vertices[ib].x() - vertices[m].x());
        }
      }
    }
  }

  return ;
}
  

template <unsigned int ORDER, typename T>
inline
const Point<SPATIAL_DIMENSION, T>* LagrangeSegmentFiniteElement<ORDER, T>::get_vertices()
{
  if (!compute_vertices_called)
    compute_vertices(vertices) ;

  return vertices ;
}


/** virtual function which evaluates the values of all
      shape functions at a given point ip */
template <unsigned int ORDER, typename T>
template <typename U>
inline
U LagrangeSegmentFiniteElement<ORDER, T>::evaluate_basis (  unsigned int ib, const Point<SPATIAL_DIMENSION, U>& px)
{
  femlib_assert(ib < Traits::N_Nodes, "Invalid basis index request");

  const Point<SPATIAL_DIMENSION, T>* vertices = LagrangeSegmentFiniteElement<ORDER, T>::get_vertices() ;

  U bvalue = 1.0 ;
  for(unsigned int j=0; j < Traits::N_Nodes; ++j)
  {
    if (j != ib)
      bvalue = bvalue * (px.x() - vertices[j].x()) / (vertices[ib].x() - vertices[j].x());
  }

  return bvalue ;
}

/** virtual function which evaluates the values of all
      shape functions at all point ips */
template<unsigned int ORDER, typename T>
template <typename U>
inline
void LagrangeSegmentFiniteElement<ORDER, T>::evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, U>& ips, std::vector<U>& basis_vals )
{
  basis_vals.resize((Traits::N_Nodes)*ips.n_points, 1.0);
  
  const Point<SPATIAL_DIMENSION, T>* vertices = LagrangeSegmentFiniteElement<ORDER, T>::get_vertices() ;

  for (unsigned int k = 0, offset=0; k < ips.n_points; k++, offset+=(Traits::N_Nodes) )
  {
    const U& x = ips.x(k) ;
    for (unsigned int i = 0; i < Traits::N_Nodes; i++ )
      for (unsigned int j = 0; j < Traits::N_Nodes; j++ )
      {
        if ( j != i )
          basis_vals[i+offset] = basis_vals[i+offset] * ( x - vertices[j].x() ) / ( vertices[i].x() - vertices[j].x() );
      }
  }
  
  return ;
}



template <unsigned int ORDER, typename T>
inline
void LagrangeSegmentFiniteElement<ORDER, T>::compute_vertices(Point<SPATIAL_DIMENSION, T> vertices[Traits::N_Nodes] )
{
  T points[Traits::N_Nodes];
  
  // eGaussLobattoLegendre
  polylib::GaussLobattoLegendre<T> glob;
  glob.compute_points_and_weights (points, weights, Traits::N_Nodes) ;

  // eGaussGaussLegendre
  // polylib::GaussLegendre<T> gleg;
  // gleg.compute_points_and_weights (points, weights, Traits::N_Nodes) ;

  for (unsigned int i=0; i < Traits::N_Nodes; ++i)
    vertices[i][0] = points[i];

  compute_vertices_called = true ;
}


template <unsigned int ORDER, typename T>
inline
void LagrangeSegmentFiniteElement<ORDER, T>::reference_element_coords(std::vector<Point<SPATIAL_DIMENSION, T> >& nodes)
{
  const Point<SPATIAL_DIMENSION, T>* vertices = LagrangeSegmentFiniteElement<ORDER, T>::get_vertices() ;

  nodes.resize(Traits::N_Nodes);

  for (unsigned int i=0; i < Traits::N_Nodes; ++i)
    nodes[i][0] = vertices[i][0] ;
}


// Some useful typedefs
typedef LagrangeSegmentFiniteElement<1, double> LinearLagrange1D ;
typedef LagrangeSegmentFiniteElement<2, double> QuadraticLagrange1D ;
typedef LagrangeSegmentFiniteElement<3, double> CubicLagrange1D ;
typedef LagrangeSegmentFiniteElement<4, double> QuarticLagrange1D ;
typedef LagrangeSegmentFiniteElement<5, double> QuinticLagrange1D ;



#endif // #ifndef __lagrange_1d_hh__
